function f = bifssgamma0(p0, gamma0Lower, gamma0Upper, theta)

g = 9.8;
m = 100;
L = 2*pi;
h = (gamma0Upper-gamma0Lower)/m;
for i = 1:m
	gamma0(i) = (i-1)*h + gamma0Lower;
	if (gamma0(i) < 0)
	f(i) = (4*pi^2*sqrt(-2*gamma0(i))/(3*L^2))*[(p0*theta/100 - p0)^3] + ...
			(-2*gamma0(i))^(1.5)*(p0*theta/100 - p0) - g*p0^2;
	else 
	f(i) = (4*pi^2*sqrt(2*gamma0(i))/(3*L^2))*[(-p0)^3-...
			(p0*theta/100 - p0)^3] + ...
			(2*gamma0(i))^(1.5)*(-p0*theta/100) - g*p0^2;
	end
end
plot(gamma0, f);
xlabel('\gamma_0', 'FontSize', 14);
title({['p0 = ', num2str(p0), ', \theta = ', num2str(theta)]});
end


